In MayaGans/CoreMicro: Testing Widely Used Core Microbiome Assignment Methods

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Load packages

library(CoreMicro)
library(patchwork)
library(ggplot2)
library(tidyverse)
library(phyloseq)

Create Human OTU and Arab OTU objects from data incldued in package

human_stool <- human_stool
names(human_stool)[1] <- "X"
human_stool$X <- as.character(human_stool$X)
human_otu <- human_stool


#read in file
arabidopsis_rhizo <- arabidopsis_rhizo
names(arabidopsis_rhizo)[1] <- "X"
arab_otu <- arabidopsis_rhizo

Figure 2

Define core_venn.core_method function for creating figure 2

core_venn.core_methods <- function(combined_taxa_data,
                                   #doing this so I can remove the figures manually. 
                                     category_names = c("a",
                                                        "b",
                                                        "c",
                                                        "d"),
                                     low = "white",
                                     high = "#7E191B") {

    temp <-
      data.frame(combined_taxa_data) %>%
      dplyr::filter(.data$value == 1) %>%
      dplyr::select(.data$X, .data$name)

    temp <- lapply(split(temp, temp$name), `[[`, "X")

    suppressMessages(
      suppressWarnings(
        ggVennDiagram::ggVennDiagram(temp, category.names = category_names, set_size = 12, label_size = 10) +
          ggplot2::scale_fill_gradient(low=low,high = high) +
          scale_color_manual(values = rep("black", 4))
      )
    )

}

calculate core venn diagram

arab_otu2<-arab_otu

arab_venn <- arab_otu2 %>% core_methods() %>% 
    core_venn.core_methods(high= "#1B301E", low = "#E1E4E1") + 
  ggtitle(expression(paste("Shared Core Membership of ", italic("Arabidopsis")))) +
  theme(plot.title = element_text(size = 20, hjust = 0.5))  +
          scale_color_manual(values = rep("black", 4))  + guides(fill="none") + ggtitle("") + theme(plot.title = element_text(size = 30))


#Add abundance information for core assignments, e.g., % core vs % non-core reads.
arab_otu$reads_per_taxa <- rowSums(arab_otu[,-1])
#assign membership to core
arab_assignments_SSR <- arab_otu %>% CoreMicro::core_methods() %>% as.data.frame() %>% dplyr::filter(name == "Summation of Sequence Reads") 
#join
arab_otu_joined_SSR <- full_join(arab_otu, arab_assignments_SSR) 

arab_assignments_HC <- arab_otu %>% CoreMicro::core_methods() %>% as.data.frame() %>% dplyr::filter(name == "Hard Cutoff") 
arab_otu_joined_HC <- full_join(arab_otu, arab_assignments_HC) 

arab_assignments_PRep <- arab_otu %>% CoreMicro::core_methods() %>% as.data.frame() %>% dplyr::filter(name == "Proportion of Sequence Replicates") 
arab_otu_joined_PRep <- full_join(arab_otu, arab_assignments_PRep) 

arab_assignments_PRR <- arab_otu %>% CoreMicro::core_methods() %>% as.data.frame() %>% dplyr::filter(name == "Proportion of Sequence Reads and Replicates") 
arab_otu_joined_PRR <- full_join(arab_otu, arab_assignments_PRR) 


stacked<-rbind(arab_otu_joined_SSR, arab_otu_joined_HC, arab_otu_joined_PRep, arab_otu_joined_PRR)

levels(stacked$name)[1] <- 'Abundance'
levels(stacked$name)[2] <- 'Occupancy'
levels(stacked$name)[3] <- 'Abundance and Occupancy'
levels(stacked$name)[4] <- 'Hard Cutoffs'

arab_abund_core_plot<-stacked %>% dplyr::select(reads_per_taxa, value, name) %>% ggplot( aes( fill=as.factor(value), y = reads_per_taxa, x = name)) +  geom_bar(position="fill", stat="identity") + ylab("Proportion of total reads") + xlab("") +  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 22))  + scale_x_discrete(labels = function(x) stringr::str_wrap(x, width = 12)) + scale_fill_discrete(name ="Core \nmembership", labels=c('non-core', 'core'))  + theme_classic() + theme(axis.text.x=element_text(size=22))+ theme(axis.text.y=element_text(size=22)) + theme( axis.title.y = element_text(size = 30))  + theme(legend.text=element_text(size=22)) + theme(legend.title = element_text(size = 24))

arab_abund_core_plot

calculate core venn diagram

human_otu2 <- human_otu

human_venn <- human_otu2 %>% core_methods() %>% 
 core_venn.core_methods(low = "#EAE6f3", high = "#432976") +
  ggtitle("") +
  theme(plot.title = element_text(size = 20, hjust = 0.5)) +
          scale_color_manual(values = rep("black", 4)) + guides(fill="none")  + theme(plot.title = element_text(size = 30)) 


#Add abundance information for core assignments, e.g., % core vs % non-core reads.
human_otu$reads_per_taxa <- rowSums(human_otu[,-1])
#assign membership to core
human_assignments_SSR <- human_otu %>% CoreMicro::core_methods() %>% as.data.frame() %>% dplyr::filter(name == "Summation of Sequence Reads") 
#join
human_otu_joined_SSR <- full_join(human_otu, human_assignments_SSR) 

human_assignments_HC <- human_otu %>% CoreMicro::core_methods() %>% as.data.frame() %>% dplyr::filter(name == "Hard Cutoff") 
human_otu_joined_HC <- full_join(human_otu, human_assignments_HC) 

human_assignments_PRep <- human_otu %>% CoreMicro::core_methods() %>% as.data.frame() %>% dplyr::filter(name == "Proportion of Sequence Replicates") 
human_otu_joined_PRep <- full_join(human_otu, human_assignments_PRep) 

human_assignments_PRR <- human_otu %>% CoreMicro::core_methods() %>% as.data.frame() %>% dplyr::filter(name == "Proportion of Sequence Reads and Replicates") 
human_otu_joined_PRR <- full_join(human_otu, human_assignments_PRR) 


stackedh<-rbind(human_otu_joined_SSR, human_otu_joined_HC, human_otu_joined_PRep, human_otu_joined_PRR)

levels(stackedh$name)[1] <- 'Abundance'
levels(stackedh$name)[2] <- 'Occupancy'
levels(stackedh$name)[3] <- 'Abundance and Occupancy'
levels(stackedh$name)[4] <- 'Hard Cutoffs'

human_abund_core_plot <- stackedh %>% dplyr::select(reads_per_taxa, value, name) %>% ggplot( aes( fill=as.factor(value), y = reads_per_taxa, x = name)) +  geom_bar(position="fill", stat="identity") + ylab("Proportion of total reads") + xlab("") +  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 22))  + scale_x_discrete(labels = function(x) stringr::str_wrap(x, width = 12)) + scale_fill_discrete(name ="Core \nmembership", labels=c('non-core', 'core'))  + theme_classic() + theme(axis.text.x=element_text(size=22))+ theme(axis.text.y=element_text(size=22)) + theme( axis.title.y = element_text(size = 30))  + theme(legend.text=element_text(size=22)) + theme(legend.title = element_text(size = 24))
#%>%

Plot figure 2 by combining both panels

figure2a <- human_venn + arab_venn #+ patchwork::plot_annotation(tag_levels = "A")#+ plot_layout(guides = "collect")
figure2b <- human_abund_core_plot + arab_abund_core_plot + plot_layout(guides = "collect")
figure2b